#include <iostream>
#include <iomanip>
#include <cmath>
#include <cstdlib>

#include "plplot.h"
#include "plplotP.h"


using namespace std;


static void hls_cmap()
{
	/* HLS
	 * 0   - red
	 * 60  - yellow
	 * 100 - green
	 * 240 - blue
	 */
	PLFLT pos[2] = {0.0, 1.0};
	PLFLT h[2]   = {240, 0};
	//	PLFLT l[2]   = {0.5, 0.5};
	//	PLFLT s[2]   = {0.9, 0.9};
	PLFLT l[2]   = {0.3, 0.6};
	PLFLT s[2]   = {1.0, 1.0};
	//	plscmap1n(256);
	plscmap1l(0, 2, pos, h, l, s, NULL);
}


void gray_cmap()
{
	PLFLT pos[2] = {0.0, 1.0};
	PLFLT r[2]   = {0.0, 1.0};
	PLFLT g[2]   = {0.0, 1.0};
	PLFLT b[2]   = {0.0, 1.0};
	plscmap1l(1, 2, pos, r, g, b, NULL);
}


static void rainbow1_cmap()
{
	//	PLFLT pos[9] = {-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0,  1.7};  // original
	//	PLFLT pos[9] = {0.0, 0.23, 0.3,  0.38, 0.45, 0.53, 0.6,  0.68, 1.0};  // normalised
	//	PLFLT pos[9] = {0.0, 0.29, 0.34, 0.4,  0.45, 0.51, 0.56, 0.62, 1.0};  // flatter
	//	PLFLT pos[9] = {0.0, 0.17, 0.26, 0.36, 0.45, 0.55, 0.64, 0.74, 1.0};  // steeper
// 	PLFLT pos[9] = {0.0, 0.11, 0.22, 0.34, 0.45, 0.57, 0.68, 0.8,  1.0};  // steeper v2

// 	PLFLT r[9]   = {0.0, 0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0, 1.0};
// 	PLFLT g[9]   = {0.0, 0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0, 1.0};
// 	PLFLT b[9]   = {0.0, 0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0, 1.0};

	PLFLT pos[7] = {0.0, 0.22, 0.34, 0.45, 0.57, 0.68, 1.0};  // steeper v2

	PLFLT r[7]   = {0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0};
	PLFLT g[7]   = {0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0};
	PLFLT b[7]   = {0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0};

	plscmap1l(1, 7, pos, r, g, b, NULL);
}

static void rainbow2_cmap()
{
	PLFLT pos[7] = {0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0};

	PLFLT r[7]   = {0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0};
	PLFLT g[7]   = {0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0};
	PLFLT b[7]   = {0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0};

	plscmap1l(1, 7, pos, r, g, b, NULL);
}


static void flame_cmap()
{
	PLFLT pos[5] = {0.0, 0.2, 0.4, 0.6, 1.0};
	PLFLT r[5]   = {0.0, 0.5, 1.0, 1.0, 1.0};
	PLFLT g[5]   = {0.0, 0.0, 0.5, 1.0, 1.0};
	PLFLT b[5]   = {0.0, 0.0, 0.0, 0.3, 1.0};

	plscmap1l(1, 5, pos, r, g, b, NULL);
}

static void test1_cmap()
{
	PLFLT pos[10] = {0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0};
	PLFLT r[10]   = {0.0, 1.0, 0.0, 0.0, 0.3, 0.8,  0.3,  1.0,  1.0,  1.0};
	PLFLT g[10]   = {0.0, 0.5, 0.4, 1.0, 0.0, 0.0,  0.2,  0.7,  1.0,  1.0};
	PLFLT b[10]   = {0.0, 0.0, 0.0, 0.0, 0.4, 1.0,  0.0,  0.0,  0.95, 1.0};

	plscmap1l(1, 10, pos, r, g, b, NULL);
}


static void aips_cmap()
{
	PLFLT pos[20] = {0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5, 0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0};
	PLFLT r[20]   = {0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
	PLFLT g[20]   = {0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8, 0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0};
	PLFLT b[20]   = {0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

	plscmap1l(1, 20, pos, r, g, b, NULL);
}
















int main(int argc, char **argv)
{
    if (argc != 2) {
        printf("need a command line arg!\n");
        exit(99);
    }
    
	int dimStart = atoi(argv[1]);
	int nSteps = 256;
//	int nSteps = dimStart;

	plsdev("xwin");
	plspage(0,0,1000,800,10,100);
	plinit();

	int dim = dimStart;
	int nx = dim;
	int ny = dim;

	double xmin = 0.0;
	double xmax = 1.0;
	double ymin = 0.0;
	double ymax = 1.0;
	double fx   = 2.0;
	double fy   = 2.0;
	double pi   = 3.141592654;
	double **data;
	double *xval = new double[nx];
	double *yval = new double[ny];

	data = new double* [nx];
	for (register int x=0;x<nx;x++)
		data[x] = new double [ny];

	double xinc = xmax/(nx-1);
	for (int x=0;x<nx;x++) {
		xval[x] = xinc * x;
	}
	double yinc = ymax/(ny-1);
	for (int y=0;y<ny;y++) {
		yval[y] = yinc * y;
	}

//	double xinc = xmax/(nx);
//	for (int x=0;x<nx;x++) {
//		xval[x] = (xinc/2.0) + xinc * x;
//	}
//	double yinc = ymax/(ny);
//	for (int y=0;y<ny;y++) {
//		yval[y] = (xinc/2.0) + yinc * y;
//	}

	double zmax = 0.0;
	double zmin = 1e10;
	for (int y=0;y<ny;y++) {
		for (int x=0;x<nx;x++) {
			//data[x][y] = x;
			//data[x][y] = xval[x];
			//data[x][y] = yval[y] + x/10.0;
			//data[x][y] = xval[x] + yval[y];
			//data[x][y] = 0.5*sin(2*pi*fx*xval[x]) + 0.5*sin(2*pi*fy*yval[y]);
			data[x][y] = 0.5*cos(2*pi*fx*xval[x]) + 0.5*cos(2*pi*fy*yval[y]);
			//data[x][y] = 10.0*cos(2*pi*fx*xval[x]);
			//data[x][y] = 2.2*cos(2*pi*fx*xval[x])*cos(2*pi*fy*yval[y]);
			if (data[x][y] > zmax) zmax = data[x][y];
			if (data[x][y] < zmin) zmin = data[x][y];
		}
	}
	cout << endl;
	cout << "===================================================" << endl;
	cout << "Data" << endl;
	cout << "===================================================" << endl;
	cout << "Data size: " << nx << " x " << ny << endl;

//	cout << endl;
//	cout << "(x, y)\t\tamp" << endl;
//	for (int y=0;y<ny;y++) {
//		for (int x=0;x<nx;x++) {
//			//cout.precision(2);
//			cout << fixed << "(" << xval[x] << ", " << yval[y] << ")\t";
//			//cout.precision(4);
//			cout << data[x][y] << endl;
//		}
//		cout << endl;
//	}
//	cout << endl;

	cout << "xinc = " << xinc << endl;
	cout << "yinc = " << yinc << endl;
	cout << "xmin = " << xmin << ", xmax = " << xmax << endl;
	cout << "ymin = " << ymin << ", ymax = " << ymax << endl;
	cout << "zmin = " << zmin << ", zmax = " << zmax << endl;


	/* Set the colour map */
//	hls_cmap();
//	gray_cmap();
//	rainbow1_cmap();
	rainbow2_cmap();
//	flame_cmap();
//	test1_cmap();
//	aips_cmap();

	pladv(0);
	plvpor(0.15,0.75,0.2,0.8);
	plwind(xmin,xmax,ymin,ymax);
	char title[50];
	plcol0(15);
	sprintf(title,"Data (%d x %d)",nx,ny);
	plmtex("t", 4.0, 0.5, 0.5, title);
	plmtex("b", 4.0, 0.5, 0.5, "x");
	plmtex("l", 4.0, 0.5, 0.5, "y");

	plimage(data,
			(PLINT) nx,   (PLINT) ny,
			(PLFLT) xmin, (PLFLT) xmax,
			(PLFLT) ymin, (PLFLT) ymax,
			(PLFLT) zmin, (PLFLT) zmax,
			(PLFLT) xmin, (PLFLT) xmax,
			(PLFLT) ymin, (PLFLT) ymax);

	plcol0(15);
	plschr (0, 0.6);
	plbox("bcnst",0.0,0,"bcnst",0.0,0);
	//plbox("bcnst",xinc,4,"bcnst",yinc,4);
	plschr (0, 1.0);







	/*======================================================================
	 * Scale
	 *======================================================================*/
	cout << endl;
	cout << "===================================================" << endl;
	cout << "Scale" << endl;
	cout << "===================================================" << endl;
	double **scale;
	int width = 2;
	scale = new double* [width];
	for (register int x=0;x<width;x++)
		scale[x] = new double [nSteps];

	double scaleInc = (zmax-zmin)/(nSteps-1);
	double scaleMin = zmin - (scaleInc/2.0);
	double scaleMax = zmax + (scaleInc/2.0);

	cout << "- zmin      = " << zmin     << endl;
	cout << "- zmax      = " << zmax     << endl;
	cout << "- scale min = " << scaleMin << endl;
	cout << "- scale max = " << scaleMax << endl;
	cout << "- scale inc = " << scaleInc << endl;

	for (int step=0;step<nSteps;step++) {
		for (int w=0;w<width;w++) {
			scale[w][step] = zmin + scaleInc * step;
		}
	}

//	cout << endl;
//	cout << "Scale: " << endl;
//	for (int step=0;step<nSteps;step++) {
//		cout << "(s" << step << ") " << scale[0][step] << endl;
//	}
//	cout << endl;

	/* Plot scale */
	pladv(1);
	plvpor(0.79,0.82,0.2,0.8);
	plwind(0.0,xmax,zmin,zmax);
	plimage(scale, 				// 2d data array
			width,nSteps,		// data dimensions 	-
			xmin,xmax,			// x plot coordinates
//			scaleMin,scaleMax,	// y plot coordinates
			zmin,zmax,
			zmin,zmax,			// z range to which to plot data
			0.0,xmax,			// x dimensions of sub window to plot
			zmin,zmax);			// y dimensions of sub window to plot
//			scaleMin,scaleMax);	// y dimensions of sub window to plot
	plcol0(15);
	plschr (0, 0.6);
	plbox("bc",0.0,0,"cmst",((zmax-zmin)/5.0),5);
//	plbox("bc",0.0,0,"cmst",0.0,0);
//	plbox("bc",0.0,0,"cmst",scaleInc,1);
	plschr (0, 1.0);
	plcol0(15);
	plschr (0, 0.5);
	plmtex("l", 2.0, 1.0, 1, "amplitude scale");
	plschr (0, 1);



	/* Free memory */
	for (register int x=nx-1; x>=0;x--)
		delete [] data[x];
	delete[] data;

	for (register int x=0; x>=0;x--)
		delete [] scale[0];
	delete[] scale;
	plend();

	return EXIT_SUCCESS;
}
